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= -T3 Method for Modeling Complex Occlusions in Fluid Simulations 

FIELD OF THE INVENTION 

This invention pertains to computational fluid dynamics, particularly for 
the purpose of generating realistic fluid animation for computer graphics. 

5 BACKGROUND 

Introduction 

An occlusion is the means by which solid objects are represented 
within a fluid simulation. A solid object is a 3D volume that is impermeable to fluid 
but whose velocity and surface attributes will affect the adjacent fluid's flow via the 
10 integration of occlusion information into the fluid simulation. 

This fluid simulation is applicable to voxel based fluid simulations 
based on finite difference methods or on finite volume methods employing the 
Navier-Stokes fluid simulation equations. 

General Fluid Simulation 
15 Eulerian-based fluid simulation consist at its core of two alternating 

steps: (1) calculate a fluid advection velocity field, which represents how the fluid 
should flow during a discrete time step, based on the current fluid properties and any 
external non-fluid influences (i.e. occlusions, gravity) via Eulerian-equations known 
to govern the behavior of fluid (i.e. the Navier-Stokes equation) and (2) advance the 
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non-velocity aspects of fluid representation across the current discrete time step 
using the newly derived fluid advection velocity field. 

Commonly, most Eulerian-based fluid flow governing equations are 
based on the Navier-Stokes equations (see equation 1). 

u t = vV • ( Vu) - (u • V)u--Vp + g 
5 P 

Equation 1 

The Navier-Stokes equations used for fluid simulation can vary but in 
general they consist of a pressure term, a convection term, a viscosity term and an 

1 o external forces term . 

The pressure projection method is the most common approach to 
solving the incompressible Navier-Stokes fluid flow equations in the context of 
computational fluid dynamics. The pressure projection method involves calculating 
first the solution to all the components of the Navier-Stokes equation except the 

15 pressure projection term and then the effect of the pressure term is integrated as a 
separate step yielding the full Navier-Stokes. This procedure is outlined in detail in 
[Stam 1999] and [Foster & Fedkiw 2001.] 

The influences of occlusions are captured in two of the above listed 
terms of the Navier-Stokes equations: the external forces term and the pressure 

20 term. The external forces term captures the contact velocity of the occlusions in 
regards to the fluid. The pressure term captures the impermeable character of 
occlusions. In the pressure projection step the impermeability of occlusions is 
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captured via the specification of von Neumann conditions on contact boundaries 
between occlusions and the fluid. 

Existing Methods 

The Figures 2, 3, 4 and 5 in this section are reproduced from [Versteeg 
5 & Malalasekera 1995.] 

Cartesian Grid Representation 

Static Occlusions 

A rough approximation to the occlusion region is achieved by denoting 
whole voxels as being either an occlusion or not. The resulting approximation is 
10 unable to represent any occlusion surface orientation that does not align with a voxel 
face. Figure 2 below illustrates the standard why in which two quarter spheres 
would be approximated on a Cartesian grid. 

Dynamic Occlusions 

For dynamic occlusions it is necessary to update the list of voxels 
1 5 denoted as occlusions at each time step over the course of the simulation. 

Integration into Fluid Simulation 

Despite their shortcomings, blocked out occlusions on Cartesian grids 
have the advantage of being quite easy to integrate into the standard Navier-Stokes 
fluid simulation equations. The velocities within occlusion cells are set to the velocity 
20 of the occlusion prior to both the advection and the pressure calculations. 
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Simulation Quality 

For occlusions whose surfaces did not originally align with any voxel 
faces the Cartesian grid approximation of the occlusion will be very poor This is 
evident in figure 3. 

5 Curvilinear Grids 

Representation 

Static Occlusions 

Curvilinear grids allow for arbitrary warping of the otherwise regular 
grid structure. Thus when occlusions are defined by denoting their component 
10 voxels the underlying grid can be warped such that voxel faces can be forced to 
align with the surfaces of the original occlusions. (See figure 4.) 

Dynamic Occlusions 

In order to ensure that voxel faces align with the surfaces of dynamic 
occlusions it is necessary to update the warping of the curvilinear grid at each time 
15 step over the course of the simulation. Because there are many fluid properties 
represented on the curvilinear grid structure it is necessary to transfer this 
information from the previous warping to the new warping as well. 

Integration Into Fluid Simulation 

Because of the arbitrary warping of the grid structure many 
20 calculations on curvilinear grids are significantly more difficult and time consuming 
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than on Cartesian grids. However, because of the significant improvement in 
occlusion representation provided by curvilinear grids they have nonetheless 
become popular. 

The quality in which occlusions are integrated into the fluid simulation 
5 when using curvilinear grids (figure 5) is significantly improved as compared to the 
quality produced from existing methods for Cartesian grids. 
SUMMARY 

According to one aspect of the present invention there is provided in a 
method of fluid simulation where state of a fluid comprising of velocities is updated in 
10 the presence of impermeable objects having surfaces in a given region over discrete 
time steps by: 

dividing the region into cells comprising a regular grid and then 
defining a velocity field which associated a velocity vector with each cell; and 

recalculating the velocity field at each consecutive time step based on 
15 the state of the fluid on the previous time step and the effect of impermeable object 
surfaces via Navier-Stokes equations comprising calculation of advection and 
pressure effects; 

the improvement comprising: 

assigning a value to the velocity vectors, associated with the cells 
20 contained within the impermeable objects when the velocity field is used for the 
calculation of the advection and pressure effects, which is copied from the closest « 
fluid containing cell; and 
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when the value includes a normal component which would cause 
motion of the fluid into the object, removing the normal component. 

When the impermeable objects have velocities defined on their 
surfaces, the method may include: 
5 determining the relative velocity by taking the difference between the 

velocity from the closest fluid containing cell and the velocity from the nearest 
impermeable object surface; and 

removing the normal component which would cause motion of the fluid 
into the object by taking the dot product of the relative velocity with the surface 
10 normal of the nearest impermeable object surface, and when it is negative, adding to 
the velocity a vector which has a magnitude of the dot product times the magnitude 
of the velocity and points in the direction of the surface normal of the nearest 
impermeable object surface. 

When a fluid volume including a surface defined by level set values 
15 representing the distance to the surface is advected according to the velocities, the 
method may include: 

storing velocity data only for those cells which are inside or near the 
fluid volume; and 

storing level set values only for those cells which are near the fluid 

20 surface. 

The method may further include obtaining the impermeable object 
surface velocities, the impermeable object surface normals; and determining 
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whether a cell is inside or outside of the impermeable objects using the level set and 
velocity field as described according to the second aspect of the present invention. 

The method may further include obtaining the velocity from the closest 
fluid cell by extrapolating the velocities from the cells just outside the impermeable 
5 object surface into the cells inside the impermeable object surface satisfying the 
constraint that the gradient of the extrapolated velocities along the direction of the 
impermeable object surface is zero. 

When the impermeable objects may be deforming and includes a 
transformation along a path, the method may include computing the impermeable 
10 object surface velocities as the sum of the velocity caused by the transformation 
along the path and the velocity caused by the deforming of the object surface. 

According to a second aspect of the present invention there is provided 
in a method of fluid simulation where state of a fluid is updated in the presence of 
impermeable objects having surfaces in a given region over discrete time steps by; 
15 dividing the region into cells comprising a regular grid; and 

recalculating the fluid state at each consecutive time step based on the 
state of the fluid on the previous time step and the effect of impermeable object 
surfaces via Navier-Stokes equations; 

the improvement comprising: 
20 defining the impermeable objects as a level set with level set values 

representing the signed distance to the nearest impermeable objects surface, in 
conjunction with a velocity field comprising the velocities of the neareast 
impermeable object surface. 
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The method may further include: 

storing level set values only for those cells which are near the 
impermeable object surface; and 

storing velocity values only for those cells which are near the 
5 impermeable object surface. 

Illustrated Example 

The following is an illustrated example of how the method described in 
this patent results in a noticeable improvement in the accuracy at which occlusions 
are represented in the fluid simulation. 

10 The test scene used is illustrated in figure 6. This scene consists of a 

long rectangular tube, open on both ends, and partially obstructed by two semi- 
elliptical occlusions. The fluid is being forced into the tube on the left open end and 
out of the right thus setting up a cross flow though the tube. 

Figure 7 illustrates the resulting fluid velocities using the standard 

15 Cartesian grid method of approximating the influence of occlusion in a fluid flow. All 
voxels containing any portion of an occlusion have been denoted (i.e. the filled block 
voxels.) Additionally, the velocities in the denoted regions have been set to the 
velocity of the containing occlusions, which in this case is zero. In this figure the 
fluid velocities resulting form this poor approximation are illustrated as the blue 

20 vectors originating from each cell center. The critical artifact is the visible 
dampening of the fluid velocities immediately adjacent to the cells denoted as 
occlusions. 
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The method described in this filing does not block out an over 
estimation of the occlusion as described above but rather under estimates the 
occlusion slightly. The velocities of occlusions cells are not set to their respective 
object velocities but rather to the velocity of the nearest fluid cell constrained such 
5 that any object-inward motion is prevented. 

It is convenient to calculate the resulting constrained velocity as the 
sum of two orthogonal components - one parallel to the occlusion normal and the 
other perpendicular: 



[proj^u, n o U/ >n u n c 
lproj B# u„ otherwise 



u ± =ortho n< (u / (l-2or u )+u o a u ) 



10 where ur is the velocities of the nearest fluid cell, u 0 Is the true velocity of the 
occlusion, a 0 is the occlusion slip coefficient, and n 0 is the surface normal of the 
nearest occlusions surface. The slip coefficient parameterizes the range of 
boundary conditions from full-slip to no-slip along the range [0.1]. 

This described method of preparing the occlusion velocities is used 

15 prior to the calculation of the advection and the pressure effects. 

The fluid velocities resulting from using this method of accurately 
incorporating the occlusions into the fluid simulation is evident in figure 8. There is 
no visible dampening of the fluid velocities directly adjacent to occlusions in this 
figure (8) unlike in figure 7. The fluid velocities in figure 8 also follow closely the 

20 smooth curved surface of the occlusions. 
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In Figure 9 is similar to figure 8 except that the black occlusion regions 

have been removed thus exposing the red underlying velocities that have been set 

using the improved description above. 

BRIEF DESCRIPTION OF THE FIGURES 
5 In the accompanying drawings, which illustrate an exemplary 

embodiment of the present invention: 

Figure 1 is an overview of the steps of the invention process: 

Figure 2 illustrates delineation of occlusions boundaries on a Cartesian 

grid. 

10 Figure 3 illustrates fluid velocities resulting from standard Cartesian 

grid occlusion representations. 

Figure 4 illustrates a curvilinear grid that incorporations occlusions; 
Figure 5 illustrates fluid velocities resulting from standard curvilinear 
grid occlusion representations. 
15 Figure 6 details the initial scene using in the example as specified on a 

Cartesian grid. 

Figure 7 illustrates fluid velocities resulting from standard Cartesian 
grid occlusion representations. 

Figure 8 illustrates fluid velocities resulting from the improved method 
20 Cartesian grid occlusion representations detailed in this document. 

Figure 9 illustrates how the occlusion detail is represented via the 
velocities interior to the occlusions (red vectors) on the Cartesian grid. 
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DETAILED DESCRIPTION 

Unless defined otherwise, all technical and scientific terms used herein 
have the same meaning as commonly understood by one of ordinary skill in the art 
to which the invention belongs. Although any methods similar or equivalent to those 
described herein can be used in the practice of the present invention, the preferred 
methods are described herein. An exemplary embodiment will now be described. 

Unified Occlusion Representation 

In order to unify all of the occlusions affecting the simulation, we 
represent their geometry using a level set instead of relying on a polygon mesh 
directly (as was done in (Foster and Fedkiw 2001) and (Enright, et al 2002)). A level 
set alone can only capture the instantaneous geometry of the occlusions, thus it 
cannot model the effects of the occlusion velocities and slip conditions. We store 
this information in two addition fields - a vector field for the velocities and a scalar 
field for the slip conditions. We denote the level set, occlusion velocity field, and slip 
condition field by the following set: 



The global occlusions representation is the union of all the individual 
occlusions. For our purposes the union operator is defined as: 



^(x)^(x),u(4<*M} 
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This union operator, since it is both associative and distributive, allows 
considerable flexibility in the calculation and partial caching of the global occlusions 
model. 

Creating Occlusions Representations 
5 The three fields that represent occlusions in our approach often need 

to be generated from explicit geometry representation. The following method is 
used to compute the set of fields for one occlusion. 

Level Set 

Of the three fields, the generation of the level set is the most complex. 

10 Our method, which is both efficient and robust, consists of the following steps: (1) 
Explicit geometry is converted into a structure that allows for efficient calculation of 
ray intersections. (2) Rays are cast along a number of axis directions, recording data 
indicating whether the center of a voxel is inside or outside, and an initial distance 
from the surface. This data is then resolved into a rough distance function. (3) 

15 Discontinuities near the zero-level set (which can result from imperfections such as 
gaps in the geometry) are removed. (4) The level set is then "smoothed" via a 
standard re-initialization algorithm. The generation of the level set is described in 
detail in prior US provisional patent application No. 60/443,543, the subject matter of 
which is incorporated herein by reference. 
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Velocity Field 

Our calculation of the velocity field is straightforward. Each occlusion 
is associated with a motion path expressible as a function mapping time t to a 
transformation matrix, T(f). The velocity field is obtained by taking the product of the 
5 position and the component-wise derivative of this motion path: u(x,f)= (T'(f))x. 

Slip Field 

We typically use a single value for the slip condition of each occlusion, 
so the slip field is a constant. It is also possible to define a varying slip field for the 
slip condition of each occlusion. 

10 Pseudocode 

ComputeUnifiedOcclusionsO { 
Set the global occlusions level set globjs to be empty. 
Set the global occlusions velocity field glob_vel to be zero. 
Set the global slip condition field glob_slip to a default value. 
15 For(each occlusion occ in the simulation) { 
Compute a level set Is from occ 
Compute a velocity field vel from occ 
Compute a slip condition field slip from occ 
For(each point x in the grid of Is) { 
20 lf(ls[x] < glob_ls(x]) { 
glob_ls[x] = ls[x] 
glob_vel[x] - vel[x] 
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glob_slip[x] = slip[x] 

} 

} 

} 

5 } 

Incorporation into Fluid Simulation 

The occlusions are relevant to two major aspects of the fluid solver, the 
mass conservation calculation and the advection velocity field. 

. Mass Conservation 

10 Incompressibility is enforced by solving the equation v .u-v' P for 

pressure, p. During the calculation of ?■„, vector occlusion velocities are used 
instead of fluid velocities in regions interior to the occlusions. The occlusions are 
incorporated into the Poisson matrix represented the Laplacian operator by omitting 
cells that lie completely within the occlusions, and using von Neumann boundary 

15 conditions for the remaining faces whose centers' lie within the occlusions. 

Advection Velocity Field 

The advection velocity field, which is used to advect both the level set 
and the fluid velocity field, must accurately account for the occlusions. This is 
achieved via our new method of constrained velocity extrapolation. This is similar to 
20 the method described in (Enright 2002) for extrapolating fluid velocities into the air, 
the main distinction however, besides the fact that we are extrapolating Into 
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occlusions rather than the air, is that constraints are imposed on the velocities after 
extrapolation. These constraints reflect the effects of three space varying aspects of 
the occlusions: the slip coefficient, the surface orientation, and the velocity. It is 
convenient to calculate the resulting constrained velocity as the sum of two 
5 orthogonal components - one parallel to the occlusion normal and the other 
perpendicular: 

u = fP«>j..«/ -.-u^-.u. Ui=orth0! ( U/ (,_2„)*„. a ) 

1 1 P r °j., °. otherwise 

Where u f is the unconstrained extrapolated fluid velocity, u 0 is the 
10 occlusion velocity, a 0 is the occlusion slip coefficient, and n 0 is the normalized 
gradient of the occlusions level set. The slip coefficient parameterizes the range of 
boundary conditions from full-slip to no-slip along the range [0,1]. 
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While one embodiment of the present invention has been described in 
the foregoing, it is to be understood that other embodiments are possible within the 
5 scope of the invention. The invention is to be considered limited solely by the scope 
of the appended Claims. 



